Quantum Monte Carlo results for bipolaron stability in quantum dots 
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Bipolaron formation in a two-dimensional lattice with harmonic confinement, representing a simplified model 
for a quantum dot, is investigated by means of quantum Monte Carlo simulations. This method treats all interac- 
tions exactly and takes into account quantum lattice fluctuations. Calculations of the bipolaron binding energy 
reveal that confinement opposes bipolaron formation for weak electron-phonon coupling, but abets a bound state 
at intermediate to strong coupling. Tuning the system from weak to strong confinement gives rise to a small 
reduction of the minimum Frohlich coupling parameter for the existence of a bound state. 
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I. INTRODUCTION 

Continuous experimental improvements in the preparation 
and investigation of semiconductor quantum dots in recent 
years have sparked a lot of interest in extending our theoretical 
understanding of charge carriers confined in such quasi-zero- 
dimensional systems. 1 Apart from potential technological ap- 
plications, quantum dots with tunable properties represent a 
playground to compare theory and experiment. 

The importance of the lattice degrees of freedom in, e.g., 
polar semiconductors has led to a large number of studies 
on polarons and bipolarons. 2 These quasiparticles correspond 
to bound states of a single carrier in a self-induced lattice 
distortion or of two electrons in a shared virtual phonon 
cloud. Bipolarons manifest themselves in, e.g., tunneling 
experiments 3 or transport through molecular dots4^ 

Whereas polaron formation in quantum dots is by now 
quite well understood (see Refs. 6,7,8,9,10 and refer- 
ences therein), conflicting results exist on the stability of 
bipolarons . 3,8 - 1Q ' 11 ' 12 i 13 - 14 These calculations are based on vari- 
ational treatments, with several works employing strong- 
coupling or adiabatic approximations. It is known from stud- 
ies of polaron and bipolaron formation that such methods are 
not able to fully capture the relevant physicS j 2 i 15 i 16 and their 
use hence represents a possible source of the contradictory 
findings. The aim of this paper is to resolve the above is- 
sues by applying an unbiased quantum Monte Carlo (QMC) 
method^ 

Previous works apply continuum models and the effective 
mass approximation. Here we argue that a lattice model is es- 
sential for several reasons. First, bipolaron physics involves 
self-trapping of carriers due to strong interaction with the lat- 
tice, a process governed by lattice fluctuations on the scale of 
the unit cell. 17 Second, for intermediate to strong coupling, 
bipolarons are rather small, so that any continuum description 
is expected to break down. Third, quantum dots studied ex- 
perimentally often contain only a relatively small number of 
unit cells. 

Bipolaron formation in a quantum dot model with local in- 
teractions has recently been studied numerically. 10 The pur- 
pose of this work is to extend these calculations to a more 
realistic model with long-range interactions — similar to the 
continuum models employed by other authors — and to obtain 
results for the bipolaron binding energy. 



II. MODEL 

The Hamiltonian considered here takes the form 

h = -t c Ls> + ^ u ^K^j,i + K J2 i ri i 2 " 1 
+ y^Z (* 2 +^ 2 ) ~ a J2h> ifiti J '' (1) 

i i,j 

with the long-range Coulomb interaction 

u ={ U ' ri = Tj (2) 

and the interaction between an electron at r, and the oscillator 
at r jt 

ki= (1^-^ + 1)3/2- ( 3 ) 

In Eq. (fl}, cj CT creates an electron with spin a at lattice site 
i (located at r^), Xi (pi) denotes the displacement (momen- 
tum) of the harmonic oscillator at site i, hi t(J — c\ a c i a 
and hi — J2 a ^i,<r- The model parameters are the nearest- 
neighbor hopping integral i, the Coulomb repulsion U, the 
confinement strength K, the (dispersionless) optical phonon 
frequency ojo and the electron-phonon coupling constant a. 
We analyze a two-dimensional (2D) square lattice with two 
electrons of opposite spin and periodic boundary conditions. 

For \r\ > 1, the Coulomb interaction Uij in Eq. (0 is of 
the same form as in previous studies of bipolarons in quantum 
dotsi 8 ' 11 ' 12 ' 14 An important difference between previous work 
and our lattice model is that does not diverge for = rj. 
Instead, there is a finite Hubbard U > for two electrons at 
the same site, which we believe to be more appropriate for a 
discussion of small bipolaron states. 

The electron-phonon interaction [Eq. (0)] can be regarded 
as a lattice version of the Frohlich mode h 18 ' 19 i 20 The Holstein- 
Hubbard model 10 is recovered by setting = USij and 

f jti = 5ij. We defined A = 2%/W and % = E P /?„, 
with the atomic-limit polaron binding energy of the Holstein- 
Hubbard model, Ep = a 2 /2ojq, and the free bandwidth 
W = 8t. In the sequel, all energies will be measured in units 
of t, and N denotes the lattice size in each dimension. We 
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further set the lattice constant to 1. The 2D case models a 
disk-shaped quantum dot with negligible thickness. The de- 
pendence on dimensionality has been studied for the Holstein- 
Hubbard model, 10 and is weak for the Frohlich model. 19 

The lattice Hamiltonian [Eq. (Q])] is appropriate in the study 
of self-trapping, strong electron-phonon coupling or strong 
confinement. Compared to other calculations of phonon ef- 
fects in quantum dots in the framework of the independent bo- 
son model^ we take into account the finite size of the dot and 
hence the possibility of electronic hopping. This is essential 
because (bi)polaron formation is determined by the balance 
between kinetic, lattice, and Coulomb energies. 

The harmonic (parabolic) confinement potential in Eq. (Q~|i 
is centered around site (0, 0) [the lattice extends from — (N — 
l)/2 to (N — 1) /2 in each dimension] and is usually assumed 
to be a good approximation for real quantum dots. 1 The use 
of optical bulk phonons is sensible as bipolaron formation is 
dominated by the coupling to such branches, 3 and the details 
of the phonon spectrum are of minor importance. 8 

The Hamiltonian [Eq. (Q])] contains the relevant terms to de- 
scribe bipolaron formation in quantum dots. Although more 
general cases can also be treated with the present method, our 
lattice model — more amenable to exact numerical treatments 
than the continuum version — is chosen to be as similar as pos- 
sible to the models in existing work on bipolaron formation in 
order to resolve the conflicts in available results. 

Finally, due to the simplicity of our model, we refrain from 
fitting our results to experiment by tuning parameters. Instead, 
we chose the adiabaticity ratio 7 = uio/t — 0.1 (i.e., in the ex- 
perimentally important adiabatic regime 7 <C 1) and U/t = 4 
as in previous work ; 10,16 ' 22 and vary the electron-phonon cou- 
pling and the confinement strength. Note that within the 
Frohlich model, U and a are not independent parameters. 18 



III. METHOD 

The worldline QMC method 10 ' 23 ' 24 can be extended to the 
case of long-range electron-electron interaction 24 ' 26 and an 
electron-phonon coupling of the Frohlich type^2i2£ Alterna- 
tive QMC schemes have also been used to study bipolaron 
formation, 

Apart from the controlled (and small) Trotter error, no 
approximations are made. In particular, the quantum- 
mechanical nature of the phonon degrees of freedom is fully 
taken into account, and all interactions are treated on the same 
footing. For the calculations, we have used a low temperature 
f3t = t/(kftT) = 15, a sufficiently small Trotter parameter 
At = 0.05, and a linear lattice size N — 31, ensuring small 
finite-size effects. As pointed out before^ ' 28 autocorrelations 
can be large and must therefore not be ignored. 

The bipolaron binding energy is defined as 

E B = E{2) - 2E{1) , (4) 

where E(N e ) denotes the total energy of the system with 
N e electrons. 31 Krishna et al^i introduced an additional es- 
timated Coulomb correlation term of two unbound polarons 



in a quantum dot, which generally leads to enhanced bind- 
ing energies. Another criterion for the stability of a bipo- 
laron based on the relative distance of the two electrons has 
also been suggested, 13 yielding a broader region of existence. 
Definition (01 has been used by most previous authors, and 
permits comparison to the case K = 0. Furthermore, Eq as 
defined by Eq. has a direct physical meaning, as it enters 
the Boltzmann factor that controls the average number of dou- 
bly occupied quantum dots in an ensemble of dots. 8 Note that 
we use the same exact method to calculate E(l) and E(2). 
We also measure the average electron-electron separation 




For A = 0, R gives an estimate of the size of the quantum dot, 
whereas for A > it may be regarded as a measure for the 
bipolaron radius. For the calculation of other observables we 
refer to a previous paper. 10 

IV. RESULTS 

A bipolaron is defined as a bound state of two electrons cou- 
pled to the lattice. The energy gain as compared to a system 
with two noninteracting polarons originates in the additional 
potential energy from sharing a common lattice polarization 
cloud. In the atomic limit t = 0, the latter increases quadrati- 
cally with the number of electrons per site. This leads to a ten- 
dency of the electrons to share the same region of space (i.e., 
to real-space pairing), and gives rise to a bound state for any 
A > if U = 0, or for A > A c (U) in the case U > where the 
electron-electron Coulomb repulsion has to be overcome^ 6 ^ 

Depending on the model and the parameters, the average 
distance between the particles can be larger than 1 (large bipo- 
laron), about 1 (intersite bipolaron) or less than 1 (small or on- 
site bipolaron). 16 For K = 0, the crossover from a state with 
two unbound polarons to a bipolaron with increasing A can be 
detected from observables such as electron-electron correla- 
tion functions or the "radius" R defined by Eq. © (see inset of 
Fig- IU- However, in a confined system, the two-particle wave 
function is "squeezed" even for A = 0, 10,13 and the binding 
energy is the only reliable indicator, with E^ < for a bound 
state. 

In some previous work on the continuum model, the 
strength of the confinement potential was measured in terms 
of a confinement length. Therefore, we start by analyzing the 
average electron-electron separation as a function of K. Fig- 
ure Q] shows results for three different values of A. For A = 
the distance rapidly decreases from R w N/2 at K = to 
R pa 2-5 at finite K. For intermediate coupling A = 0.5 the 
dependence is similar, but R is systematically smaller than for 
A = due to the bipolaron effect. Finally, for strong coupling 
A = 1, a small bipolaron is the ground state even for K = 0, 
so that confinement has very little influence on R. The range 
of R in Fig. Q] is comparable to the confinement lengths stud- 
ied by other authors, 14 and the different behavior for weak and 
strong couplings is consistent with continuum calculations. 8 
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FIG. 1: Average electron-electron distance R as a function of con- 
finement strength K for different values of the electron-phonon cou- 
pling parameter A. The inset shows J? as a function of A for two 
values of K. Here 7 = 0.1, U/t = 4, pt = 15 (At = 0.05), and 
N — 31. Lines are guides to the eye, and error bars are smaller than 
the symbols. 

The direct contribution of the confinement term in Eq. (Q~|i 
to the total energy is linear in N e , and therefore cancels when 
calculating Eg. Consequently, any effect of K on Eg is indi- 
rect, mediated by changes in the interaction energies related 
to lattice distortions and Coulomb correlation. 

In Fig. 13 a) we plot the bipolaron binding energy Eg as a 
function of A for different values of if. As expected, all curves 
are strictly monotonic decreasing with increasing electron- 
phonon coupling. The binding energy is positive (unbound 
state) for A < A c , and negative for A > A c (bound state or 
bipolaron). The inset in Fig. [2 a) reveals a reduction of A c with 
increasing K by about 10%, and the K = critical coupling 
is similar to that of the Holstein-Hubbard model [cf. Fig. 7(a) 
in Ref. 10]. Of course A c depends on U, with A c — > as 
U — > 0. For strong coupling, all curves in Fig.[2|a) eventually 
collapse due to the formation of a small-bipolaron state which 
is rather insensitive to confinement. 

From Fig. |2a), we can identify three different regimes, 
for which we plot the binding energy as a function of if in 
Fig. 0b). In the regime A -C A c , Coulomb interaction is 
stronger than the phonon-mediated electron-electron attrac- 
tion, so that we have two polarons with Eg > 0. Confine- 
ment enhances both interactions, and our results reveal that 
the Coulomb energy becomes even more dominant, i.e., Eg 
increases with increasing K for K/t < 1. This is particularly 
evident for the case A — shown in Fig. |2|a), where Eg in- 
creases noticeably with K, approaching the value of the onsite 
Coulomb repulsion U/t = 4. 

For A 3> A c , a bound state exists as a result of the lat- 
tice deformation energy winning over the Coulomb repulsion, 
and confinement acts in favor of bipolaron binding by further 
increasing |£b|- The influence of K on the results is much 
weaker than that for A <C A c because the physics is dominated 
by local correlations, as reflected by the small bipolaron ra- 
dius (see inset of Fig. [TJ). The weak dependence of A c on K 
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FIG. 2: (Color online) Bipolaron binding energy (a) Eb as a function 
of A for different K, and (b) Eb as a function of K for different A. 
The dashed horizontal line indicates Eb = 0, and the inset in (a) 
shows a closeup view. 



may be attributed to the fact that confinement increases both 
Coulomb repulsion as well as the phonon-mediated attractive 
interaction, and thus leaves the effective interaction and A c al- 
most unchanged. 

In between these two limits, for intermediate A, Eg shows a 
nonmonotonic behavior as a function of K. For intermediate 
coupling, we find a maximum in Eg at small values of K/t, 
and a decrease for an even stronger confinement. Whereas for 
A = 0.5 Eg > for all K shown, a bound state arises for a 
strong-enough confinement K/t > 1 and A = 0.6. In both 
cases, the variation of Eg with K is rather small. For large 
K, the curves in Fig.[2fb) approach the corresponding values 
of the effective onsite interaction U — 2Ep which plays an 
important role in molecular quantum dots with only a single 
electronic level4^ 

The main findings of previous calculations for the contin- 
uum Frohlich model are as follows. 81013 There exists a min- 
imal X C (U) below which a bound state is suppressed due to 
Coulomb repulsion. Moreover, some authors argue that there 
is also a critical confinement strength beyond which the di- 
verging Coulomb repulsion suppresses a bound state for any 
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finite electron-phonon coupling strength A < oo i ' • No 
unanimous conclusion has been reached concerning the effect 
of (weak to intermediate) confinement on the size of the pa- 
rameter region of existence for the bound bipolaron stated 

Let us relate these findings to our results. Similar to previ- 
ous work and models without confinement, we find a critical 
coupling A c for the formation of a bound state. As pointed out 
before, within our lattice model, A c changes slightly as a func- 
tion of the confinement strength K, and scales approximately 
linearly with l/M 

Concerning the regime of strong confinement, our choice 
of a finite Coulomb repulsion of two electrons with opposite 
spin located at the same lattice site leads to substantially dif- 
ferent physics. Even for K — > oo, a bound onsite bipolaron 
is formed if the effective onsite interaction U — 2Ep < 0. Fi- 
nally, we observe that confinement opposes binding for weak 
coupling A <C A c , whereas it enhances binding for intermedi- 
ate to strong coupling A > A c . Despite the associated changes 
of the value of the binding energy, the phase diagram (i.e., the 
region of existence of a bound state) is only weakly affected 
by the confinement potential. This insensitivity is even more 
pronounced for the Frohlich parameter a c oc \/A^, which de- 
creases with increasing K/t € [0, 2] by about 5%. 

Previous authors have used material-specific parameters to 
make predictions for the existence of bipolarons in typical 
quantum dot systems. 14 As we believe that our model is too 
simple to make quantitative statements, we restrict ourselves 
to a mostly qualitative discussion. The (small) reduction of the 
critical coupling due to confinement suggests that the dot size 
may determine whether bipolarons are stable or not. However, 
according to our findings, this should only be relevant for ma- 
terials with intermediate electron-phonon coupling. In con- 
trast, for weak coupling, bipolaron formation is suppressed by 
Coulomb repulsion, whereas for strong coupling a bipolaron 



ground state is stable regardless of the confinement strength. 
Besides, as demonstrated in Fig. |2{b), the value of the bind- 
ing energy does have a noticeable dependence on the confine- 
ment potential, at least for experimentally relevant weak to 
intermediate values of the electron-phonon coupling. The dot 
size may also change the phonon spectrum or electronic band 
structure, but such effects have been neglected in our model. 

Finally, we can estimate the order of magnitude of the bipo- 
laron binding energy in real systems by assuming a typical 
bandwidth of 1 eV. From Fig. [2] we conclude that forplausible 
values of 0.5 < A < 1, Eb is a fraction of an eV for U = 0.5 
eV, a value falling into the narrow-band regime U t. The 
binding energy increases with decreasing U, but the physics of 
bipolaron formation remains qualitatively the same. Our val- 
ues of Eg are larger than previous variational 3D results^ 4 - Ex- 
perimentally, apart from pair tunneling, 3,5 the bipolaron effect 
should also manifest itself in shot-noise measurementSi^Ml 

In summary, we have presented unbiased quantum Monte 
Carlo results for bipolaron formation in a two-dimensional 
quantum dot, taking into account the crystal lattice and quan- 
tum phonon effects. Confinement is found to give rise 
to — respectively strengthen — a bound state at intermediate to 
strong electron-phonon interaction, and to reduce the critical 
coupling for bipolaron formation. The present method can be 
used to study more general models with dispersive phonons or 
more complicated electronic bands and dot geometriesJSi^ 5 . 
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